home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / intpol.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  3.8 KB  |  115 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. ;;; Interpolation routine by CFFK.
  13. (macsyma-module intpol)
  14. (load-macsyma-macros transm numerm)
  15.  
  16. (declare-top (special $intpolrel $intpolabs $intpolerror)
  17.      (flonum $intpolrel $intpolabs a b c fa fb fc)
  18.      (fixnum lin)
  19.      (notype (interpolate-check flonum flonum flonum flonum))) 
  20.  
  21. (COMMENT  |  For historical information ONLY.  |
  22. (defun fmeval2 (x) 
  23.        (cond ((integerp (setq x (meval x))) (float x))
  24.          ((floatp x) x)
  25.          (t (displa x) (MAXIMA-ERROR '|not floating point|))))
  26. (defun qeval (y x z) (cond (x (fmeval2 (list '($ev) y (list '(mequal) x z) '$numer)))
  27.                (t (funcall y z))))
  28. )
  29.  
  30. (or (boundp '$intpolabs) (setq $intpolabs 0.0)) 
  31. (or (boundp '$intpolrel) (setq $intpolrel 0.0))
  32. (or (boundp '$intpolerror) (setq $intpolerror t))
  33.  
  34. (Defun $interpolate_SUBR (F LEFT RIGHT)
  35.   (BIND-TRAMP1$
  36.    F F
  37.    (prog (a b c fa fb fc (lin 0))
  38.      (declare (flonum a b c fa fb fc) (fixnum lin))
  39.      (setq A (FLOAT LEFT)
  40.        B (FLOAT RIGHT))
  41.      (or (> b a) (setq a (prog2 niL b (setq b a))))
  42.      (setq fa (FCALL$ f a)
  43.        fb (FCALL$ f b))
  44.      (or (> (abs fa) $intpolabs) (return a))
  45.      (or (> (abs fb) $intpolabs) (return b))
  46.      (and (> (*$ fa fb) 0.0)
  47.       (cond ((eq $intpolerror t)
  48.          (merror "function has same sign at endpoints~%~M"
  49.              `((mlist)
  50.                ((mequal) ((f) ,a) ,fa)
  51.                ((mequal) ((f) ,b) ,fb))))
  52.         (t (return $intpolerror))))
  53.      (and (> fa 0.0)
  54.       (setq fa (prog2 nil fb (setq fb fa)) a (prog2 nil b (setq b a))))
  55.      (setq lin 0.)
  56.     binary
  57.      (setq c (//$ (+$ a b) 2.0)
  58.        fc
  59.        (FCALL$ f c))
  60.      (and (interpolate-check a c b fc) (return c))
  61.      (cond ((< (abs (-$ fc (//$ (+$ fa fb) 2.0))) (*$ 0.1 (-$ fb fa)))
  62.         (setq lin (f1+ lin)))
  63.        (t (setq lin 0.)))
  64.      (cond ((> fc 0.0) (setq fb fc b c)) (t (setq fa fc a c)))
  65.      (or (= lin 3.) (go binary))
  66.     falsi
  67.      (setq c (cond ((> (+$ fb fa) 0.0)
  68.             (+$ a (*$ (-$ b a) (//$ fa (-$ fa fb)))))
  69.            (t (+$ b (*$ (-$ a b) (//$ fb (-$ fb fa)))))) 
  70.        fc (FCALL$ f c))
  71.      (and (interpolate-check a c b fc) (return c))
  72.      (cond ((> fc 0.0) (setq fb fc b c)) (t (setq fa fc a c)))
  73.      (go falsi))))
  74.  
  75. (defun interpolate-check (a c b fc)
  76.        (not (and (prog2 nil (> (abs fc) $intpolabs) (setq fc (max (abs a) (abs b))))
  77.          (> (abs (-$ b c)) (*$ $intpolrel fc))
  78.          (> (abs (-$ c a)) (*$ $intpolrel fc)))))
  79.  
  80.  
  81.  
  82.  
  83. (DEFUN INTERPOLATE-MACRO (FORM TRANSLP)
  84.        (SETQ FORM (CDR FORM))
  85.        (COND ((= (LENGTH FORM) 3)
  86.           (COND (TRANSLP
  87.              `(($INTERPOLATE_SUBR) ,@FORM))
  88.             (T
  89.              `((MPROG) ((MLIST) ((msetq) $NUMER T))
  90.                    (($INTERPOLATE_SUBR)  ,@FORM)))))
  91.          ((= (LENGTH FORM) 4)
  92.           (LET (((EXP VAR . BNDS) FORM))
  93.            (SETQ EXP (SUB ($LHS EXP) ($RHS EXP)))
  94.            (COND (TRANSLP
  95.               `(($INTERPOLATE_SUBR)
  96.                 ((LAMBDA-I) ((MLIST) ,VAR)
  97.                     (($MODEDECLARE) ,VAR $FLOAT)
  98.                     ,EXP)
  99.                 ,@BNDS))
  100.              (T
  101.               `((MPROG) ((MLIST) ((msetq) $NUMER T))
  102.                     (($INTERPOLATE_SUBR)
  103.                      ((LAMBDA) ((MLIST) ,VAR) ,EXP)
  104.                      ,@BNDS))))))
  105.          (T (merror "wrong number of args to INTERPOLATE"))))
  106.  
  107. (DEFMSPEC $INTERPOLATE (FORM)
  108.   (MEVAL (INTERPOLATE-MACRO FORM NIL)))
  109.  
  110. (def-translate-property $INTERPOLATE (FORM)
  111.   (let (($tr_numer t))
  112.     (TRANSLATE (INTERPOLATE-MACRO FORM t))))
  113.  
  114.  
  115.